fsave <- function(x, file, location = "./data/processed/", ...) {
if (!dir.exists(location))
dir.create(location)
datename <- substr(gsub("[:-]", "", Sys.time()), 1, 8)
totalname <- paste(location, datename, file, sep = "")
print(paste("SAVED: ", totalname, sep = ""))
save(x, file = totalname)
}
fpackage.check <- function(packages) {
lapply(packages, FUN = function(x) {
if (!require(x, character.only = TRUE)) {
install.packages(x, dependencies = TRUE)
library(x, character.only = TRUE)
}
})
}
colorize <- function(x, color) {
sprintf("<span style='color: %s;'>%s</span>", color, x)
}
## Loading required package: ggmap
## Google's Terms of Service: https://cloud.google.com/maps-platform/terms/.
## Please cite ggmap if you use it! See citation("ggmap") for details.
## Loading required package: leaflet
## Loading required package: tmap
## Loading required package: basemapR
## [[1]]
## NULL
##
## [[2]]
## NULL
##
## [[3]]
## NULL
##
## [[4]]
## NULL
##
## [[5]]
## NULL
##
## [[6]]
## NULL
# Loading 100m-by-100m raster data:
rast <- sf::st_read(dsn = "./data/rawGIS/cbs_vk100_2021_v1.gpkg")
## Reading layer `vierkant_100m_2021' from data source
## `/Users/robsalasco/Proyectos/bigssslabjournal/data/rawGIS/cbs_vk100_2021_v1.gpkg'
## using driver `GPKG'
## Simple feature collection with 386024 features and 34 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 13900 ymin: 307500 xmax: 277500 ymax: 611500
## Projected CRS: Amersfoort / RD New
# Next we load the shapefile of the administrative neighbourhoods ('buurt') and districts ('wijk'):
neighbShape <- sf::st_read(dsn = "./data/rawGIS", layer = "buurt_2021_v1")
## Reading layer `buurt_2021_v1' from data source
## `/Users/robsalasco/Proyectos/bigssslabjournal/data/rawGIS' using driver `ESRI Shapefile'
## Simple feature collection with 14175 features and 43 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 10425.16 ymin: 306846.2 xmax: 278026.1 ymax: 621876.3
## Projected CRS: Amersfoort / RD New
districtShape <- sf::st_read(dsn = "./data/rawGIS", layer = "wijk_2021_v1")
## Reading layer `wijk_2021_v1' from data source
## `/Users/robsalasco/Proyectos/bigssslabjournal/data/rawGIS' using driver `ESRI Shapefile'
## Simple feature collection with 3331 features and 40 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 10425.16 ymin: 306846.2 xmax: 278026.1 ymax: 621876.3
## Projected CRS: Amersfoort / RD New
# ... And then the zipcode shapes:
postcode4Shape <- sf::st_read(dsn = "./data/rawGIS", layer = "CBS_pc4_2020_v1")
## Reading layer `CBS_pc4_2020_v1' from data source
## `/Users/robsalasco/Proyectos/bigssslabjournal/data/rawGIS' using driver `ESRI Shapefile'
## Simple feature collection with 4068 features and 36 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 13565.4 ymin: 306846.2 xmax: 278026.1 ymax: 615045.3
## Projected CRS: Amersfoort / RD New
postcode5Shape <- sf::st_read(dsn = "./data/rawGIS", layer = "CBS_pc5_2020_v1")
## Reading layer `CBS_pc5_2020_v1' from data source
## `/Users/robsalasco/Proyectos/bigssslabjournal/data/rawGIS' using driver `ESRI Shapefile'
## Simple feature collection with 33287 features and 36 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 13565.4 ymin: 306846.2 xmax: 278026.1 ymax: 615045.3
## Projected CRS: Amersfoort / RD New
postcode6Shape <- sf::st_read(dsn = "./data/rawGIS", layer = "CBS_pc6_2020_v1")
## Reading layer `CBS_pc6_2020_v1' from data source
## `/Users/robsalasco/Proyectos/bigssslabjournal/data/rawGIS' using driver `ESRI Shapefile'
## Simple feature collection with 460392 features and 34 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 13565.4 ymin: 306846.2 xmax: 278026.1 ymax: 615045.3
## Projected CRS: Amersfoort / RD New
# If we want to remove them:
rast <- rast[rast$aantal_inwoners != -99997, ]
# If we want to replace -99997 with an arbitrary integer in [1,4]:
# rast$aantal_inwoners[rast$aantal_inwoners == -99997] <- 2
# If we want to replace -99997 with NA: rast$aantal_inwoners[rast$aantal_inwoners == -99997] <- NA
rast <- sf::st_transform(x = rast, crs = sf::st_crs("+proj=longlat +datum=WGS84"))
rast <- sf::st_centroid(rast)
## Warning in st_centroid.sf(rast): st_centroid assumes attributes are constant over geometries of x
neighbShape <- sf::st_transform(x = neighbShape, crs = sf::st_crs("+proj=longlat +datum=WGS84"))
districtShape <- sf::st_transform(x = districtShape, crs = sf::st_crs("+proj=longlat +datum=WGS84"))
postcode4Shape <- sf::st_transform(x = postcode4Shape, crs = sf::st_crs("+proj=longlat +datum=WGS84"))
postcode5Shape <- sf::st_transform(x = postcode5Shape, crs = sf::st_crs("+proj=longlat +datum=WGS84"))
postcode6Shape <- sf::st_transform(x = postcode6Shape, crs = sf::st_crs("+proj=longlat +datum=WGS84"))
city <- "Gouda"
shape <- districtShape[districtShape$GM_NAAM == city,]
shape$color <- sample(rainbow(n = nrow(shape), alpha = 1))
# Specifying a "blank" theme for our ggplot visualisation:
mapTheme <- ggplot2::theme(
legend.position = "bottom",
panel.border = element_blank(),
axis.title = element_blank(),
axis.ticks = element_blank(),
axis.text = element_blank()
)
plot <- ggplot() +
base_map(st_bbox(shape), increase_zoom = 2, basemap = 'positron') +
geom_sf(data = shape, aes(fill = color), alpha = .7) +
guides(fill = FALSE) +
mapTheme
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
plot

# Overlaying the district shapes:
#plot <- ggmap(mapTiles, darken = c(0.8, "white")) +
#ggtitle(city) +
# geom_sf(
# data = shape,
# aes(fill = WK_NAAM),
# color = "black",
# alpha = 0.5,
# inherit.aes = FALSE
# ) +
# coord_sf(datum = NA)
#print(plot + mapTheme + theme(legend.position = "none"))
#
labels <- sf::st_centroid(shape) |> sf::st_coordinates() |> as.data.frame()
## Warning in st_centroid.sf(shape): st_centroid assumes attributes are constant over geometries of x
labels$label <- shape$WK_NAAM
plotlabels <- geom_text(
data = labels,
#position = position_dodge2(0.1), # option to reduce overplotting
aes(x = X, y = Y, label = label)
)
print(plot + plotlabels + mapTheme + theme(legend.position = "none"))

city <- "Rotterdam"
shape <- subset(
districtShape,
districtShape$GM_NAAM == city & districtShape$AANT_INW >= 0 # removes NA's
)
tm_shape(shape) +
tm_fill(col="AANT_INW", title = "Population", id="WK_NAAM", style="jenks")

zoom = 12 # Higher values give more detailed basemaps. 12-15 usually suffice.
coords <- as.data.frame(sf::st_coordinates(shape))
mapTiles <- ggmap::get_stamenmap(
bbox = c(
left = min(coords$X),
bottom = min(coords$Y),
right = max(coords$X),
top = max(coords$Y)),
maptype = "toner-background",
crop = FALSE, zoom = zoom
)
## Source : http://tile.stamen.com/toner-background/12/2093/1352.png
## Source : http://tile.stamen.com/toner-background/12/2094/1352.png
## Source : http://tile.stamen.com/toner-background/12/2095/1352.png
## Source : http://tile.stamen.com/toner-background/12/2096/1352.png
## Source : http://tile.stamen.com/toner-background/12/2097/1352.png
## Source : http://tile.stamen.com/toner-background/12/2098/1352.png
## Source : http://tile.stamen.com/toner-background/12/2099/1352.png
## Source : http://tile.stamen.com/toner-background/12/2100/1352.png
## Source : http://tile.stamen.com/toner-background/12/2093/1353.png
## Source : http://tile.stamen.com/toner-background/12/2094/1353.png
## Source : http://tile.stamen.com/toner-background/12/2095/1353.png
## Source : http://tile.stamen.com/toner-background/12/2096/1353.png
## Source : http://tile.stamen.com/toner-background/12/2097/1353.png
## Source : http://tile.stamen.com/toner-background/12/2098/1353.png
## Source : http://tile.stamen.com/toner-background/12/2099/1353.png
## Source : http://tile.stamen.com/toner-background/12/2100/1353.png
## Source : http://tile.stamen.com/toner-background/12/2093/1354.png
## Source : http://tile.stamen.com/toner-background/12/2094/1354.png
## Source : http://tile.stamen.com/toner-background/12/2095/1354.png
## Source : http://tile.stamen.com/toner-background/12/2096/1354.png
## Source : http://tile.stamen.com/toner-background/12/2097/1354.png
## Source : http://tile.stamen.com/toner-background/12/2098/1354.png
## Source : http://tile.stamen.com/toner-background/12/2099/1354.png
## Source : http://tile.stamen.com/toner-background/12/2100/1354.png
## Source : http://tile.stamen.com/toner-background/12/2093/1355.png
## Source : http://tile.stamen.com/toner-background/12/2094/1355.png
## Source : http://tile.stamen.com/toner-background/12/2095/1355.png
## Source : http://tile.stamen.com/toner-background/12/2096/1355.png
## Source : http://tile.stamen.com/toner-background/12/2097/1355.png
## Source : http://tile.stamen.com/toner-background/12/2098/1355.png
## Source : http://tile.stamen.com/toner-background/12/2099/1355.png
## Source : http://tile.stamen.com/toner-background/12/2100/1355.png
mapTheme <- ggplot2::theme(
legend.position = "left",
panel.border = element_blank(),
axis.title = element_blank(),
axis.ticks = element_blank(),
axis.text = element_blank()
)
ggmap(mapTiles, darken = c(0.8, "white")) +
ggtitle(city) +
geom_sf(
data = shape,
aes(fill = AANT_INW),#, label = WK_NAAM),
color = NA,
alpha = 0.5,
inherit.aes = FALSE
) +
scale_fill_viridis_c() +
labs(fill = "Population") +
coord_sf(datum = NA) +
mapTheme
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.

suppressMessages(sf_use_s2(FALSE))
# Adding administrative area information:
rast <- sf::st_intersection(x = sf::st_as_sf(rast), y = neighbShape)[, c(1:39, 78)] # selecting only relevant columns
## although coordinates are longitude/latitude, st_intersection assumes that they are planar
## Warning: attribute variables are assumed to be spatially constant throughout all geometries
# Adding postcode information:
rast <- sf::st_intersection(x = sf::st_as_sf(rast), y = postcode6Shape)[, c(1:40,74)] # selecting only relevant columns
## although coordinates are longitude/latitude, st_intersection assumes that they are planar
## Warning: attribute variables are assumed to be spatially constant throughout all geometries
# We now have the 6-digits postcodes; the 4- and 5- digits postcodes then are:
rast$PC5 <- substr(rast$PC6, start = 1, stop = 5)
rast$PC4 <- substr(rast$PC6, start = 1, stop = 4)
fsave(rast, compress = TRUE, file = "raster.RData")
## [1] "SAVED: ./data/processed/20220727raster.RData"
# Loading the file
load("./data/processed/20220712raster.RData") # Make sure filename is correct!
rast <- x
rm(x)
# Plotting the raster of a city:
city = "Rotterdam"
cityrast <- rast[rast$GM_NAAM == city,]
palette <- leaflet::colorFactor(
palette = "Set3", # or other RColorBrewer palettes e.g "Greens", "magma"
domain = cityrast$BU_NAAM
)
leaflet::leaflet(cityrast) |>
leaflet::addTiles() |>
leaflet::addProviderTiles(providers$Stamen.Toner) |>
leaflet::addCircles(
label = ~BU_NAAM,
color = ~palette(BU_NAAM),
opacity = 0.7
)
## Warning in RColorBrewer::brewer.pal(max(3, n), palette): n too large, allowed maximum for palette Set3 is 12
## Returning the palette you asked for with that many colors
## Warning in RColorBrewer::brewer.pal(max(3, n), palette): n too large, allowed maximum for palette Set3 is 12
## Returning the palette you asked for with that many colors
load("./data/processed/20220708polling_df")
pollstations <- x
rm(x)
load("./data/processed/20220712raster.RData")
rast <- x
rm(x)
pollstations <- sf::st_geometry(sf::st_as_sf(x = as.data.frame(pollstations), crs = sf::st_crs("+proj=longlat +datum=WGS84"),
coords = c("long", "lat")))
# head(pollstations) #just some points, the coordinates of the polling stations
voronoi <- sf::st_voronoi(x = pollstations %>% st_combine(), envelope = NULL)
## Warning in st_voronoi.sfc(x = pollstations %>% st_combine(), envelope = NULL): st_voronoi does not
## correctly triangulate longitude/latitude data
#This ensures that 'voronoi' has the correct CRS
voronoi <- sf::st_sf(
sf::st_collection_extract(voronoi, type = "POLYGON"),
crs = sf::st_crs("+proj=longlat +datum=WGS84")
)
#This will be the "id" of each Voronoi tile:
voronoi$voronoi <- 1:nrow(voronoi)
leaflet::leaflet(voronoi) |>
leaflet::addTiles() |>
leaflet::addProviderTiles(providers$Stamen.Toner) |>
leaflet::addPolygons(color = "blue") |>
leaflet::addCircles(data = pollstations, color = "red") |>
leaflet::setView( # This defaults the map view so that it points to Amsterdam
lng = 4.9041,
lat = 52.3676,
zoom = 14
)
suppressMessages(sf_use_s2(FALSE))
rast <- sf::st_intersection(x = rast, y = voronoi)
## although coordinates are longitude/latitude, st_intersection assumes that they are planar
## Warning: attribute variables are assumed to be spatially constant throughout all geometries
city = "Rotterdam"
cityrast <- rast[rast$GM_NAAM == city, ]
palette <- leaflet::colorFactor(sample(colors(), length(unique(cityrast$voronoi))), domain = cityrast$voronoi)
leaflet::leaflet(voronoi) |>
leaflet::addTiles() |>
leaflet::addProviderTiles(providers$Stamen.Toner) |>
leaflet::addPolygons(color = "blue") |>
leaflet::addCircles(data = pollstations, color = "red", opacity = 1) |>
leaflet::addCircles(data = cityrast, label = ~voronoi, color = ~palette(as.factor(voronoi)), opacity = 0.7) |>
leaflet::setView(lng = 4.9041, lat = 52.3676, zoom = 14)
fsave(rast, compress = TRUE, file = "raster_vor.RData")
## [1] "SAVED: ./data/processed/20220727raster_vor.RData"
neighbShape <- neighbShape[, c(1, 2, 42:44)]
districtShape <- districtShape[, c(1:4, 39:41)]
postcode4Shape <- postcode4Shape[, c(1, 37)]
postcode5Shape <- postcode5Shape[, c(1, 37)]
postcode6Shape <- postcode6Shape[, c(1, 35)]
fsave(list(neighbShape, districtShape, postcode4Shape, postcode5Shape, postcode6Shape, voronoi), compress = TRUE,
file = "shapes.RData")
## [1] "SAVED: ./data/processed/20220727shapes.RData"